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Abstract 

A numerical technique that solves the parabolized form of the Navier-Stokes equa- 
tions is presented. Such a method makes it possible to obtain very detailed descriptions 
of the flowfield in a relatively modest CPU time. The present approach is based on 
a space-marching technique, uses a finite volume discretization and an upwind flux- 
difference splitting scheme for the evaluation of the inviscid fluxes. Second order ac- 
curacy is achieved following the guidelines of the the ENO schemes. The methodology 
is used to investigate three-dimensional supersonic viscous flows over symmetric cor- 
ners. Primary and secondary streamwise vortical structures embedded in the boundary 
layer and originated by the interaction with shock waves are detected and studied. For 
purpose of validation, results are compared with experimental data extracted from 
literature. The agreement is found to be satisfactory. In conclusion, the numerical 
method proposed seems to be promising as it permits, at a reasonable computational 
expense, investigation of complex three-dimensional flowfields in great detail. 
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1 Introduction 


In the study of three-dimensional supersonic flows it is common to be faced with complex 
interactions concerning shock waves and viscous layers. Such occurrences often provoke 
dramatic changes in the flowfield features, both qualitatively and quantitatively. Shock- 
induced separations of the viscous layer make the flow vortex-dominated close to the wall, 
and vortical structures are also likely to appear, for sufficiently high Reynolds numbers, 
in zones where convective effects are preponderant, due the instability of the slip surfaces 
resulting from shock/shock interactions. Important consequences of these interactions are 
increases in heat fluxes, skin friction coefficients and pressures at the wall in correspondence 
with the reattachment of the separated flow; in addition transition, shock wave shapes and 
the efficiency of the air-intakes that might swallow such streams are affected. 

Numerical tests are usually helpful in the investigations of such complicated fluid-dynamic 
patterns, as they provide a tool to observe and possibly to understand the origins and 
the effects of the numerous phenomena triggered by shock/shock and shock/viscous layer 
interactions. It is clear, however, that in order to obtain good numerical results comparable 
with experimental data, a fairly detailed description of the flowfield is necessary, especially 
when multiple vortical structures are present. 

The only completely correct way of solving numerically three-dimensional compressible 
viscous flows is to integrate in time the full Navier-Stokes equations until a steady-state (if 
one exists) is reached. This approach is certainly affordable today, but, if many grid points 
are needed to solve in detail complex fluid-dynamic features, it could be excessively time 
and memory consuming. In the case of supersonic steady-state flows, however, this practical 
difficulty can be partially circumvented with the aid of the approximate form of the full 
Navier-Stokes equations known as Parabolized Navier-Stokes equations. 

As will be pointed out in the next sections, the advantage of the Parabolized Navier- 
Stokes (PNS) equations is that they can be solved using a space-marching technique, a 
characteristic which allows one to spend relatively short computational effort and also results 
in noticeable memory savings. Therefore, it is possible to reinvest time and memory in 
more refined grids, thus permitting a better resolution of the flowfield. As a drawback, 
the parabolizing assumption requires the freestream Mach number to be supersonic and the 
streamwise velocity to be always positive (streamwise flow separations are thus excluded, 
while crossflow separations are permitted); moreover, the streamwise pressure gradient must 
be altered in the subsonic part of the flowfield [1], 

In the approach presented here, the governing equations are integrated in an explicit fash- 
ion and the physical domain is discretized according to a finite volume technique. The con- 
vective part of the equations (inviscid fluxes) is treated following a flux-difference-splitting 
method with an approximate solution of a Riemann problem at each cell interface [2] [3] 
while the diffusive terms (viscous fluxes) are calculated using a centered scheme. Second 
order accuracy is achieved by means of an Essentially Non Oscillatory scheme [4] with linear 
reconstruction of the solution at each step of integration. Presently, only inert gases in lam- 
inar regime are considered, but a future extension to include thermochemical or turbulence 
effects is certainly possible. 

To validate the method, numerical results are compared with experimental data ex- 
tracted from the literature. Supersonic corner flows configurations have been chosen as a 
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benchmark as they give the chance to have, on the same geometry, shock/shock interactions 
and shock/viscous layer interactions with multiple separations. The comparison appears to 
be satisfactory, as the same flowfield features are recognized, and measured and computed 
values show a good agreement. 

2 Governing equations 

2.1 Starting point: the three-dimensional Navier-Stokes equations 

Compressible viscous flows are governed by the Navier-Stokes equations, that written in 
integral conservative form read like: 

J- / W dV + f Fi ■ ndS + / F v • ndS = 0 (1) 

at Jv Js Js 

where V represents an arbitrary volume inclosed in a surface S with unit normal ti positive 
if directed outward. 
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Figure 1: Control volume 


System (1) can be reduced to the non-dimension al form with the help of the following 
reference values: L for length, for temperature, y/RTZ for velocity, RT^ for energy per 
unit mass and //(Too) for viscosity. Therefore, from now on, the flowfield variables should be 
considered as non-dimensional. In particular, W is the vector of conservative variables 

W = {t>,pq,E} T 

tensor Fj contains the inviscid fluxes 

Fi = {pq, p7 + pq® q, ( E + p)q] T 


and tensor Fy contains the viscous fluxes 


Fv = ,-f,-kVT-T-q} T 

Quantities p , p and q = {u, v , w} T are respectively the local density, pressure and velocity; 
E represents the total energy per unit volume: 

' i«n 


E = p ( e + 


( 2 ) 
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e is the internal energy per unit mass, and Re^ are the freestream Mach number and 
Reynolds number, 7 is the ratio of the specific heats and, finally, / is the unit matrix. 
Viscous stresses are contained in tensor r, with 


Tij = ft 


'dqi.dqA 2 

,97, + aZ -3 (v<r)Sii 


(3) 


where 6,j is the Kroneker’s symbol. The viscosity is calculated via Sutherland’s law: 

'l +T ref \ 


where 


rji 3/2 


T re j = 


K T + T re f 

110.4/v 


and the thermal conductivity k is obtained according to the relation: 

H 7 


k = 


Pr 7 — 1 


(4) 


(5) 


(6) 


where Pr is the Prandtl number. Finally, the perfect gas relationship completes the set of 
equations 

71 

(7) 


P - = T 


2.2 Three-dimensional Parabolized Navier-Stokes equations 

The three-dimensional Parabolized Navier-Stokes equations are derived from the steady-state 
full Navier-Stokes equations with the aim of obtaining a system of equations representing a 
well posed problem with respect to an integration performed using a space-marching tech- 
nique. For this purpose, any derivative in the streamwise direction contained in the stress 
tensor is neglected, all viscous and heat fluxes in the streamwise direction are dropped and 
the pressure gradient in the subsonic layer is properly altered. With such modifications, 
which are valid only for sufficiently high Reynolds numbers, system (1) is reduced to a set 
of hyperbolic-parabolic equations [1][5]. 

Thus, assuming the streamwise direction to be coincident with the x-axis, the components 
of the stress tensor necessary to evaluate the viscous fluxes are reduced to the following form: 


2 , 

(8) 

2 /„ 

-n {2v y - w z ) 

(9) 

2 

— // (2 W z Vy) 

(10) 

T yx = VU y 

(11) 

T zx = 

(12) 

T zy = (t{Vz + W y ) 

(13) 
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Moreover, the x-component of the temperature gradient is neglected: 

vr = T y j + T z k ( 14 ) 

The streamwise pressure gradient is split according to the technique suggested by Vi- 
gneron, Rakich and Tannehill in reference [1]: 

Px = up x + (1 - w)p* (15) 


If only the first term of the RHS of equation (15) is retained and the second term is considered 
as a source term, the set of equations containing only the convective terms is hyperbolic 
provided that 


u > 0 


U < 


7 Mj 

i + (7 - mi 


for M . r < 1 


(16) 


where M x is the local Mach number in the streamwise direction. The first condition contained 
in equations (16) prevents the PNS approximation from being used when the flow separates 
in the streamwise direction. A short analysis of the second condition reveals that u> is equal 
to 1 when M x = 1 and is zero when M x = 0; this means that the effect of the streamwise 
pressure gradient is completely neglected at the wall, but is progressively included as the 
flow approaches supersonic conditions. For M x greater than 1 the effect of p x will be entirely 
taken into account and the value of w will be unity. Of course, since a space-marching 
integration is desired, the freestream Mach number will necessarily have to be supersonic. 

Combining the cited assumptions, system (1) is finally reduced to the following form: 


where 


and /* and /** 


/ FJ • ndS + j F* v n*dS = j s P • ndS + P' 
F| = {pq,pT+ pq®q,[E + p)q} T 

FV = 

-tc^QC 

P = {0, — (1 — cu) pi**, 0} T 
are scalar matrices with: 


(17) 


(18) 


diag I* = (u;, 1, 1) 
diag /** = (1,0,0) 

Since the viscous and heat fluxes in the steamwise direction are dropped, vector n* contains 
only the components in the y-direction and in the z-direction of the normal unit vector, that 
is n* = (0,n y ,n 2 ). 

A source term P' must be added to make the integral formulation coherent with equation 
(15), as noticed in references [6] and [7]. In fact, the corresponding integral form of up x is: 

oj J pn x dS = J <jopn x dS — p J u>n x dS (19) 
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0 z 


Figure 2: Computational control volume 

Therefore, we have: 

P' = {0,0, 0} T (20) 

where 0 is a scalar matrix with 

diag 0 == (p j urixdS, 0, 0^ (21) 

Real gas effects, that can arise for Mach numbers above 6, are neglected here. Moreover, 
the flowfield regime is assumed to be laminar. 

3 Numerical solution of the PNS equations 

3.1 Discretization and integration 

Being a set of hyperbolic-parabolic equations, PNS equations can be integrated in the stream- 
wise direction according to a space-marching technique. In the present approach, such a 
direction is assumed to be coincident with the x-axis, so that the flowfield is solved step by 
step in planes normal to the x-direction. The physical domain is discretized according to 
a finite volume approximation. Computational control volumes, as shown in figure 2, are 
hexahedrons having two faces (labeled 2 and 4) normal to the x-direction. Keeping in mind 
the previous discussion, system (17) can be discretized in the following way: 

t + t ( F v)i ' S, - •£ (n = 0 (22) 

i= 1 i= 1 «=1 

where i is the surface label. 

In equation (22), (Ff) 4 is the solution obtained at the previous integration step and 

[(Fi), - (P' ) 2 j is the unknown. The remaining inviscid and viscous fluxes across the lateral 
surfaces have to be estimated and this will be the subject of the next two subsections. 
The source term P has been omitted in equation (22), this implies that only a part of the 
streamwise pressure gradient is considered in the subsonic layer. The reason for neglecting 
P is that it should be evaluated using a backward difference (i.e. considering the previous 
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volume) and this, as explained in reference [1], would arise in a stability condition imposing 
a lower bound for Ax. 

The integration of equation (22) is performed in an explicit manner, so that the values of 
the primitives variables that are used to evaluate the lateral fluxes are those that have been 
already calculated at the previous step, that is on surface 4. The x-momentum component 
of P', which is the only one different from zero (see equations (20) and (21)), is discretized 
as: 


j P 2 = P4 (w4n* 4 AS4 + uqrixj A<Si + u) 3 n X3 AS 3 + u}$n Xi A«5s + u)qti x& £±Sq + A^) — 

= (p') 4 + tn + (p% + (p') 5 + (p\ + (p')= 

The unknown fluxes [(FJ) 2 — (P') 2 ] are finally evaluated from equation (22): 


[(FJ), - (P') 2 j • n 2 A5 2 =- [(F;)„ ■ n t ASt - (P') 4 + 

(Fj), • n,A 5, + (FV), • n*i AS, - (P'), + 
(FJ) 3 ■ n 3 AS 3 + (FV) 3 • n-jAS. - (P') 3 + 
(F) ) 5 • n s A«S s + (Fy) s ■ n‘sAS$ — (P') 5 + 
(FJ) 6 • n 6 AS 6 + (FV) 6 • n'oASe - (P') 6 ’ 


(23) 


To terminate the step, it is necessary to extract the primitive variables from the solution 
fluxes [(FJ) 2 — (P') 2 i which contain u. To accomplish this, a cubic function u — u> (F| — P') 
has been obtained starting from the cubic function u> = uj (Fj) suggested by J. Korte in [6]. 
Denoting the jth component of [(FJ) 2 — (P') 2 | f° r brevity by ( Fj) a , the form of the function 


is: 


(1 + Bp4)u>l — 


2(1 + 07 


7 


1 


+ 2 B 


7^P4 

,7-1 


- (F2)s 


U>2 + 


(1 + 0 7 2 + A _ 4 


with 


(7-ir 




7-1' 


U?2 — 


2x4^7 

(7 -l) 3 


= 0 


(24) 


7 2 (F 2 ) 


2(F 5 )JF 1 ) s -(F 3 )1-(F 4 ) 


B = 


7 2 P4 


(-, - l) 2 [2 (F s ), (Fr), - (F 3 )= - (F 4 )J] 


and where p 4 is the pressure on surface 4. The quantity c is a safety factor (~ 0.8) which is 
applied to w, that is evaluated at first according to equation (16), but is then corrected as: 

u — min (l,cw) 

Equation (24) is solved using a Newton-Raphson iteration. Once u> 2 is known, it is possible 
to evaluate Ff 2 and finally the primitive variables on the face labeled with 2. 
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3.2 Evaluation of the lateral convective fluxes 

Lateral inviscid fluxes are evaluated defining and solving an appropriate Riemann problem 
across each lateral surface (1,3,5 or 6 ). 

The definition of the Riemann problem consists, at first, in considering the two finite 
volumes connected by the lateral surface and in fixing a direction 77 joining the volumes and 
belonging to the y-z plane that contains surfaces labeled as 4. In the present approach 77 
is normal to the trace of the considered surface in the y-z plane (see figure 3). Then, the 
variation of the flowfield variables along 77 is to be considered. Due to the discretization, 
two piecewise constant (first order accuracy) or piecewise linear (second order accuracy) 
distributions of the flowfield variables are present between cells A and B, separated by a 
discontinuity in correspondence of the lateral surface (see figures 3 and 4). The collapse of 
such a discontinuity generates a pattern of waves along which signals propagate. The waves 
split the domain in the vicinity of the discontinuity in a set of uniform regions where the 
values of the flowfield variables are to be found, generating in this way a Riemann problem. 

To obtain waves directions and corresponding signals, the equations governing the invis- 
cid part of the flowfield are written in quasi-linear form in a new local frame of reference 
constituted by direction 77 , by the x-direction and by a (-direction normal to plane 77 — x (see 
figure 3). Here an approximate solution of the Riemann problem is sought for [2], and conse- 
quently the Euler equations are written under the assumption of isentropic flow. Moreover, 
the Riemann problem is solved for simplicity in two dimensions rather than in three, so that 
only variations of the flowfield variables along the x -77 plane are considered. Therefore, all 
the derivatives in the (-direction are neglected, and finally the following set of equations is 
obtained: 

r kP x + aP n + = 0 

a 2 

<7 X (X&n "I" n ( P rj du. 1 P x ) — 0 

7 u z 

< W x -f d"Wr, = 0 (25) 

2 

h x + ahj, - — ( uP x + aP v ) = 0 
. K + vh° = 0 



Figure 3: Local frame of reference 
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f 


N,M 


N.M+1/2 


N,M+I 


surface 


Figure 4: Piecewise constant distribution along rj between cells A and B 
where P, d, w, h , and h° have been chosen as dependent variables, with: 

K = 7 (l-u;)+u;^l-^j 

P = lnp 
d = rjycr + t) z t 
w = (yV + C z w 


h° = h + 


u 2 + V 2 + w 1 


2 

V w 

<7 = — T = — 
U U 


and where f} = ( rj y , tj z ) and ( = {( y , ( z ) are unit vectors defining directions 77 and (. 

Owing to their hyperbolic nature, the quasilinear equations can be replaced by the com- 
patibility equations that describe the convection of signals. In the present case, the collapse 
of the discontinuity between A and B generates a pattern of five waves (characteristic lines) 
along which the signals (Riemann invariants) defined by the compatibility equations propa- 
gate. The Riemann invariants and characteristic slopes corresponding to system (25) are: 


= dP + 


7 u 


^a 2 fl — u 2 d(a — 1 ) 

^ 2 = *_^ji+4V + “ V(1 ““ ) 


7 1 + ud 2 
dR 3 — dw = 0 
dRi = dh° = 0 


1 + U }& 2 


da = 0 


(26) 


and 


with 


Ai,5 — 


au 2 d ^ erfi 


^ 2 , 3)4 — d 


u (u 2 — a 2 ) + 7 u 2 (1 — u>) 


(27) 




a = ^[l+u; + 7 (l -w)] 
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surface 



Figure 5: The Riemann problem 


Equation (27) shows that there are only three distinct waves, which in the following will 
be called I (wave #1), II (waves #2,3,4) and III (wave #5). Such waves generate two new 
uniform regions c and d , which add up to the initial ones a and b (see figure 5). As signals 
traveling on each wave are not altered when crossing a characteristic of a different family, it 
is possible to write: 


dR\ b 

= dR\ d 



dR j 2 a 

— dRi c 

dR 2 b 

= dR2 d 

dR^a 

= dR 3c 

dRs b 

= dR 3d 

dR 4a 

— dR 4c 

dR 4b 

= dR 4d 

dRb a 

= dR 5c 




(28) 


Moreover, across the contact surface II the pressure and the streamline slope a remain 
constant, and therefore: 

d c = bi P c = Pd (29) 


To calculate the values of the flowfield variables in regions c and d, equations (28) must be 
integrated. Following the technique presented in reference [2], the integration is performed 
in an approximate way; therefore, the complete set of equations that it is necessary to solve 


is: 


R$ a — Rh c 

7?2 0 — 7?2 c 

7?2 6 = f?2 d 

Rib = Rid 

w a = w c 


~ P c + 
~ h c - 

~ hi - 

~ Pd - 


a 2 /3 — u 2 


« 2 *(a- l)y a c 


aw 


(1 + *) 


] H wcr 2 


Pc + 


vrc (1 — w) 
1 + wa 2 


a 2 w (1 + 

,7 1+^/5 
' 111 2 
K a 2 (3 + u 2 a (a — 1) 


-) 

>* 2 )i 


. u 2 d( 1 — w)\ _ 
Pd + — r - : — I a d 


1 +wa 2 


<?d 


(30) 


m = w d 

K = K 
K = K 

- Pc = Pd 

&c — <?c 
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Figure 6: The sonii line 

- - - w ----- 

Once regions c and d are known, it is possible to verify if the acoustic waves I and III are 
expansion or compression fans (that is shock waves)T It is also possible to verify in which 
region the trace of the lateral surface in the 77 -x planeTs embedded, and therefore to calculate 
the correct fluxes. In the situation depicted in figure_5, for example, the lateral surface is in 
region d, so that the flux across it can be easily foundjusing the values computed there. This 
also corresponds to splitting the flux difference between A and B in three contributions: 

(F«» - (Fi), = [(FI). - (F|)J + [(FJ),“- (FJ)J + [(F[) d - (F;),1 

and to considering the flux across the lateral surface as constituted of the flux in a plus the 
contribution of the left traveling waves: 1 

(FiW,/ 2 = (Pi). + [(Fi) i -(Ft)J + [(Fi) c -(Fi)J = (F;) J 

or of the flux in b plus the contribution of the right traveling waves: 


(Filv,,.,/! = (Fa + [( f;) 4 ^if;),] = (f;), 

In case the lateral surface is embedded in an expansion or compression fan, the values of 
the flowfield variables in correspondence with the ’sonic point’ labeled with a star in figure 
6 must be known. According to what has been previously stated, the flow here is obtainable 
from conditions: - 


R, ~ P* - 

6 \a 2 /3 + u 2 


R 2b ~h*- 


iv b = W 

K = h°' 

\l = \ 


«! \ y 

a 2 u; (1 + pm + ^u 2 cr(l -u) 


7 1 + LOG 2 


+ LOG 


surface 


(31) 


System (31) can be easily solved through iterations. When the values of the dependent 
variables in correspondence with the sonic line are known, it is possible to write, in case of 
an expansion fan: 

(fiwi /2 = (FD, + [(Fir - (Fi)j = (Fir 
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or 

(FJW+./2 = (*!). + [(*1)' - (Fi).l = (FJ)‘ 

and in case of a compression fan: 

(Fiw +1/2 = (Fa + [(fhj - (Fin 

or 

(FfW +1/a = (Fi)„ + [(FJ), - (11)1 + K*1)i - (F|)J = (FJ), + (FJ) d - (F|)* 

The approximation that the acoustic waves are isentropic also in case they are shocks 
could lead to some inaccuracies when strong shock waves are concerned. For this reason, 
when the presence of a strong compression is recognized through an acoustic wave, an exact 
Riemann problem solver is switched on. In the case of figure 6, this consists in imposing 
the Rankine-Hugoniot conditions across wave III, the Prandtl-Meyer conditions across wave 
I and the constancy of p and b across wave II. 

3.3 Evaluation of the lateral diffusive fluxes 

To evaluate the diffusive fluxes contained in vector Fy it is necessary to compute it, v and 
w and the gradients of u, v, w and T in the y-z plane in correspondence with the trace of 
the considered lateral surface. 

Gradients are calculated using the Gauss’ theorem: 

J V/ dV = J^f ■ ndS (32) 

and applying it to the y-z plane containing surfaces labeled with 4. The discretized form of 
equation (32) reduced to two dimensions is: 

fy ~ ^2 fi n y,A£i 

f. 

9 surface center 
X averaged values 
O streamwise derivative 

computational cell surface 

viscous cell 

0 Z 

Figure 7: The viscous cell 
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where A £,■ indicates each side of the viscous cell that has to be built up around the trace of 
the considered lateral surface. 

The values of /, are computed as the averaged value of / itself in that point. With 
reference to figure 7, we have for instance: 


(/») N,M+ 1/2 
(/.) jV,Af+l/2 


— ^ — (f 12 A£ 12y + /23A<?23j, + /34A44j, + /4lA*41,) 

^1234 X 

77 — — (fl2^12 z + /23 A^23i + fz4^M z + /ll A£jlJ 


where 

/l2 = J + //V-l,Af+l + Jn,M+ 1 + /N,Af) 

/23 = In,M + 1 

/34 = £ {f N,M + fN,M+l + //V+1,M+1 + fN+l,M) 

/41 = Jn,m 

On the other hand, the values of u, v and w in correspondence with the trace of the 
lateral surface are simply computed as the averages of the values in the adjacent cells: 


In,m+ 1/2 — 


Jn,m + In,m+i 
2 


3.4 Boundary conditions 

Boundary conditions have been implemented imposing no slip and adiabatic or isothermal 
conditions at the wall and, for what concerns the pressure gradient, reducing Navier-Stokes 
equations to the boundary layer limit: 

q w = 0 

f (VT) W — 0 for adiabatic wall 
1 T w imposed for isothermal wall 

(Vp)u, = 0 

3.5 Initial conditions 

The use of a space marching technique requires the specification of the flow conditions 
across the inlet surface. The geometrical configurations considered here posses a leading 
edge sufficiently thin to produce a shock wave attached to the body, permitting us to start 
the integration directly from the leading edge, where freestream conditions are imposed. It 
is necessary, however, to start with short integration steps, to allow the viscous layer to 
establish as quickly as possible. 
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Figure 8: The ENO scheme 


3.6 Second order accuracy 


Second order accuracy is achieved following the guidelines of the Essentially Non Oscillatory 
schemes for shock capturing techniques, originally proposed for time-dependent problems [4J. 
In the present case, the initial data distribution of the primitive variables U = {p, u, v, w, h} 
along the rj direction is considered to be piecewise linear. With reference to figure 8, it is for 
instance: 

U = Va + sa{t} — t}a) in A 

U = Ug + sb(t] — r) B ) in B 

The slopes s are chosen in the ENO fashion to avoid spurious oscillations at flow disconti- 
nuities, that is, for each element s^ of vector s^: 


$4 = minmod = 



. o 


^>0 


. e 


< 0 


(33) 


with 


4 - 


Ub-U a 

d(B , A) 
Ua-Uc 


y A — 


d(A , C) 

The discontinuity originating the Riemann problem is thus defined by the values: 


(Ua- i jw+1/2) l — U,4 + S A~^—^ — — — U>1 + S a 

(Uj^TV+l^)^ — Ub — S a — „ — — Ufi — Sa 


djA, B) 
d(A, B ) 


2 2 

Second order accuracy in the hyperbolic space marching direction x is obtained by using 
the (complete) quasilinear form of the governing equations, which makes it possible to obtain 
the derivatives p x , u Xl v x , w x and h x . The necessary derivatives in the y and z direction are 
evaluated using the Gauss’ theorem as in the computation of the viscous fluxes. Therefore, 
the final values for defining the Riemann problem are provided by: 


+ (U,),, — 


(UAf,jvr+i/2) L = U.4 + sa ^ 

(u n , m+1/! ) k = Ub - Sb 'MiEl + (Ux)b “ 
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3.7 Evaluation of the integration step size 

The amplitude of the integration step is chosen according to the inviscid CFL condition 
modified with a viscous correction. The adopted procedure is the following: 

Ax = min (Axj,,Ax 2 ) 


with 


Axj, 


= K, 


Ay 

I Imax 


A.X 2 — I\.vig 


Az 

1^2 I max 


„ _ 1 

hms - 2 y 

Re p u min (Ay, Az) 

4 Results and Conclusions 

To verify the capacity of the above discussed methodology to solve complex flowfields, a 
numerical study of symmetric supersonic corner flows in laminar regime is presented. The 
obtained results are also compared with experimental data taken from the literature, in 
particular from reference [8]. 

To match with experiments, the following geometrical and fluid-dynamic conditions have 
been considered: 
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where the meaning of 6 and 6 is shown in figure 9. 

The fluid-dynamic pattern typical of this kind of flows consists of a system of five shock 
waves. Two of them, which separate regions I and II in figure 10, are generated by the pres- 
ence of the wedges; the remaining three are due to the presence of a Mach disk (irregular 
reflection) in the region where the previous two converge. Contact surfaces directed towards 
the symmetry plane are generated at points where shocks interact (triple points) because 
of the different levels of entropy produced by the wave system on either side of the inter- 
action. Shock waves separating regions II and IV impinge on the boundary layer and are 
reflected as expansion waves, which in turn encounter the slip surfaces and are transformed 
in compression waves. The interaction between the impinging shock and the boundary layer 
provokes the separation of the latter in the crosswise direction. Therefore, a streamwise 
vortex develops, resulting as an obstacle to the crossflow, and thus generating a compression 
analogous to the one typical of two dimensional supersonic flows over a ramp. 
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In the following, numerical results obtained using a stretched computational grid com- 
posed of 100x100 cells will be shown. The stretching function used is: 


€n — Zw + (£ max ) 


1 -& + 


1 + 


2fh 

' & + 1 ’ 

A-K 


(34) 


with 


0 = 


\ A max ~ 1 ) 


In particular, a stretching parameter f3 y = 0 Z = 1.03 has been used. 

The external boundaries have been defined starting with y max — 0.055 and z max = 0.055 
at the inlet surface; proceeding with the integration, the upper and lateral boundary surfaces 
have been kept parallel to the opposite walls. 

The presence of the vortical structures and of the related compression fan transports the 
three-dimensional effects due to the corner at noticeable distances from it inside the shock 
layer, making it necessary to fix rather ample external boundaries. If a square or rectangular 
grid is used, this results in a lot of points outside the shock layer. To avoid wasting time 
in computing points for which it is known a priori that the freestream conditions apply, a 
procedure has been implemented that detects the position of the shock layer at each step of 
integration and solves the flowfield only in the part of the domain containing it. 

In figure 11a, crossflow streamlines corresponding to x=0.09 m are plotted. The presence 
of a large streamwise vortex split in two parts that are tied by a saddle point can be clearly 
seen. Below it, close to the wall, a small secondary vortex is captured thanks to the great 
refinement of the mesh. This picture is enlarged in figure 11c. In figure lib, pressure 
contours are presented using the same scale. The patterns described by the sketch of figure 
10 are found. 

In figure 12, limiting streamlines at the wall are plotted; crossflow streamlines have been 
placed beside for a better comprehension of the picture. The separation and reattachment of 
the primary double vortex are well evidenced, and also the presence of the secondary vortex 
at the wall can be perceived. 

The same computed results presented above are now compared with experimental data 
extracted from reference [8]. In figure 13, pitot pressure contours corresponding to a distance 
from the leading edge of 0.09 m are shown. The contours show a qualitative good agreement. 
In fact, the wedge-shock position is not exactly the same, as the computed one corresponds 
to a local slope of about 4.9° with respect to the wedge, versus an angle of about 5.5° for the 
experiments (the inviscid value should be 3.3°). Nevertheless, the maximum pitot pressure 
value is 6.70 for the computation, and experiments show their highest value on the contour 
corresponding to 6.42. 

In addition, in figure 14a, numerical and experimental static pressure at the wall are 
overlayed, showing a satisfactory agreement. In figure 14b, the flow direction at the wall is 
compared. The two curves are very similar close to the corner, as demonstrated by the fact 
that the reattachment of the primary vortex is detected in the same position. Nevertheless, 
away from the intersection of the ramps, some not negligible differences start to appear. In 
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the computation, the secondary vortex is smaller and closer to the corner (though not to such 
a great extent) and the separation of the primary vortex is located in a different position 
with respect to the experiments. Such discrepancies might possibly be due to the fact that, 
approaching the lateral boundaries of the model, experiments are affected by side-effects 
or by interference with the tunnel boundary layer. Finally, in figure 14c, the heat flux at 
the wall is compared. As it could be expected after the previous discussion, the results are 
similar as far as the reattachment of the primary vortex is concerned, with the peak of heat 
transfer in the same position and showing a similar magnitude. Conversely, the location of 
the peak related to the secondary vortex is different, though the value is the same. Other 
discrepancies concern the fact that, in the numerical results, we find no trace of the peak 
which, in the experimental curve, is signaled at y ~ 0.09; in effect, it seems difficult to 
justify its presence, since no visible reattachment is detected in that position, neither in the 
experiment. Last, it can be noticed that the computation predicts an almost null heat flux 
locally at the corner, in contrast with the experiments; in this case we think it is reasonable 
to trust the numerical result, since very close to the corner the temperature varies very 
smoothly, and on the other hand the measurement technique used cannot approach close to 
the corner. 

In order to adequately address the accuracy of the computation, the results obtained with 
the 100x100 grid are compared with those resulting from a 60x60 grid. The main features of 
the flowfield are unchanged, as can be seen from figure 15, where pitot pressure contours are 
shown. The major difference is related to the secondary vortex: with the finer mesh it is fully 
captured, while with the coarser one its presence is just sensed. A proof to this statement 
can be found in figures 16 and 17. In figure 16, crossflow streamlines are presented: in the 
left picture (60x60 mesh), the undulation of the crossflow streamlines suggests the presence 
of the secondary vortex, but in the right one (100x100 mesh), the vortex is evident. In figure 
17, the flow direction at the wall is shown again: it can be seen that in the case of the finer 
mesh (100x100), a crossflow reversal is present at about y — 0.12, while with the coarser 
mesh (60x60) this feature does not exist. 

In conclusion, the agreement arising from the comparison appears to be satisfactory, 
though some minor discrepancies appear. 

Therefore, the numerical method proposed seems to be promising, as it permits one to 
investigate complex three-dimensional flowfields in great detail of a reasonable computational 
cost. 
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Figure 12: Limiting streamlines at the wall and crossflow streamlines. 
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Figure 17: Flow direction at the wall: (a)60x60 mesh, (b)100xl00 mesh. 
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